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The description of transport phenomena in devices consisting of arrays of tunnel junctions, and 
the experimental confirmation of these predictions is one of the great successes of mesoscopic physics. 
The aim of this paper is to give a self-consistent review of sequential transport processes in such 
devices, based on the so-called "orthodox" model. We calculate numerically the current- voltage 
(7-V^) curves, the conductance versus bias voltage (G-V) curves, and the associated thermal trans- 
port in symmetric and asymmetric two-junction arrays such as Coulomb-blockade thermometers 
(CBTs), superconducting- insulator- normal- insulator-superconducting (SINIS) structures, and su- 
perconducting single-electron transistors (SETs). We investigate the behavior of these systems at 
the singularity-matching bias points, the dependence of microrefrigeration effects on the charging 
energy of the island, and the effect of a finite superconducting gap on Coulomb-blockade thermom- 
etry. 



I. INTRODUCTION 



Quasiparticle transport processes across metallic junctions play a fundamental role in the functioning of many 
devices used nowadays in mesoscopic physics. One such device is the single electron transistor (SET), invented and 
fabricated almost two decades ago [H, which has found remarkable applications as ultrasensitive charge detector Q 
and as amplifier operating at the quantum limit Q . In the emerging field of quantum computing the superconducting 
SET has been proposed and used as a quantum bit With advancements in lithography techniques, this device can 
be fully suspended thus providing a new avenue for nanoelectromechanics. Similar devices are currently used (in 
practice arrays with several junctions turn out to provide a larger signal-to-noise ratio) as Coulomb blockade (CBT) 
primary thermometers Q. Also, superconducting double-junction systems with appropriate bias can be operated as 
microcoolers The functioning of these three classes of devices is based on the interplay between two out of the 
three relevant energy scales: the superconducting gap, the charging energy, and the temperature. 

For example, in the case of microcoolers, the temperature and the gap are finite, and the charging energy is typically 
zero. A natural question to raise is then what happens if the charging energy is no longer negligible, for example if one 
wishes to miniaturize further these devices. In contradistinction, for CBTs the charging energy and the temperature 
are important, and the superconducting gap is a nuisance. A solution is to suppress the gap by using external 
magnetic fields, an idea which makes these temperature sensors more bulky and risky to use near magnetic-field 
sensitive components. Therefore, understanding the corrections introduced by the superconducting gap could provide 
an interesting alternative route, although, with present technology, the level of control required of the gap value could 
be very difficult to achieve. Finally, electrometers and superconducting SETs are operated at low temperatures, with 
the charging energy and the gap being dominant. However, large charging energies are not always easy to obtain for 
some materials due to technological limitations, while achieving effective very low electronic temperatures is limited 
by various nonequilibrium processes. 

In this article we present a unified treatment of these three devices by solving the transport problem in the most 
general case, when all three energy scales are present. Our goal is to give an eye guidance for the experimentalist 
working in the field, showing what are the main characteristics visible in the I-Vs and G-Vs, resulting from sequential 
tunneling. Both the electrical and the thermal transport are calculated in the framework of a generalized so-called 
"orthodox" theory, which includes the superconducting gap. For completeness, we offer a self-consistent review of this 
theory, which has become nowadays the standard model for describing sequential quasiparticle transport processes 
in these devices. We ignore Josephson effects which, for the purpose of this analysis just add certain well-known 
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features at certain values of the bias (the Joscphson supercurrent at zero bias and the Josephson-quasiparticle (JQP) 
peak at 2A for symmetric SETs) to the characteristic. Higher-order tunneling phenomena such as Andreev reflections 
and cotunneling are neglected as well (the former is visible only for junctions with high transparencies, which is not 
our case, the latter is a small, second-order transport effect across the island). Also the effect of the environment 
is not discussed here, since it depends on the specifics of the sample (for example the coupling capacitances of the 
electrodes). But even stripped down to the essentials of sequential tunneling, the physics of these devices is complex 
enough to justify a careful analysis, as we attempt below. We present in detail the theory and the numerical methods 
that can be used to characterize these devices. Also, a lot of interest exists nowadays in the study of hybrid structures 
- in which the electrodes and the island are made from different materials, thus having different superconducting gaps. 
Various technologies have been developed to fabricate single-electron transistors with superconducting materials other 
than Al, most notably with Nb which has a much larger gap 0. We extend our analysis also to such hybrid SETs. 

This paper is organized in the following way: besides the introductory part, there are two main sections. Section 
im and Section IIIIl followed by closing remarks in Section IIVI Section [II] discusses tunneling of electrons from one 
superconductor into another through an insulating barrier. Here we show how to calculate the currents that flow 
through the junction by using the Fermi golden rule fSubsection IH A]) . and we give a comprehensive presentation of 
the numerical methods used to calculate the integrals over the density of states and Fermi functions which appear 
naturally in this theory (Subsection IH B]) . In Section IHII wc use these results to describe transport phenomena in 
two-junction devices. We assume a quasi-cquilibrium situation (electron distribution described by a Fermi function 
with a certain effective temperature) for both the island and the leads. We start with a presentation of the master 
equation in Subsection IHI Al A first application of this formalism is described in Subsection IHI Bl where we analyze 
the effect of a finite superconducting gap for the functioning of Coulomb-blockade thermometers. Then we discuss 
the effect of a nonzero charging energy for superconductor-insulator-normal metal-insulator-superconductor (SINIS) 
structures fSubsection IHI C) ) - which normally have a negligible charging energy. After that, in Subsection IHI Dl we 
consider the case of superconducting single-electron transistors with either identical materials in the leads and island 
(symmetric SET) or with an island with a gap different from that of the electrodes (asymmetric SET). The effect 
of the gate in such transistors is further investigated in Subsection IHI El We end this section with a calculation of 
expected cooling effects in such structures (Subsection IHI F]) . We delegate some of the calculations of integrals in 
Appendix and of expansions in Appendix [BJ 

II. TRANSPORT IN SINGLE TUNNEL JUNCTIONS 

Standard lithographic techniques enable the fabrication of junctions with almost any geometrical characteristics 
down to the linewidth of the respective system. When using scanning-electron microscopes, the typical linewidths 
obtained are of the order of 20 nm and even below. More advanced methods, such as focused ion beam lithography, 
have started to be used recently. The metals are evaporated in an ultra-high vacuum evaporator, using the shadow 
evaporation technique. This technique consists of evaporation through a mask realized in the resist, at different 
angles, so that some lines would overlap in a certain region, thus defining the area of the junction. In between 
two evaporations, the sample is oxidized (typically aluminum oxide), a process which results in the formation of the 
insulating barrier. 

In this section, we assume that the area of the junction thus formed is large enough so that capacitive effects can 
be neglected (Coulomb-blockade effects are discussed in the next section). This allows us to concentrate on tunneling 
effects only. The only energy scales relevant here will then be the temperature and the tunneling matrix element. 

A. General theory of tunneling 

Consider the case of two conductors (leads) that are connected through a thin insulating barrier. The setup is 
illustrated in Fig. [T] This creates effectively a potential barrier between the "left" and "right" electrodes: however 
electrons can still tunnel from one side to another via the tunneling effect. The total many-body Hamiltonian of the 
system comprises then two parts: one is the Hamiltonian corresponding to the left and right electrodes (which can be 
either the free-electron Hamiltonian in the case of metals, or the BCS Hamiltonian in the case of superconductors), 
the other is the tunneling Hamiltonian. This Hamiltonian was first introduced by Bardeen Q: 

Ht^ rfc„fc.4.c,^.+h-c-, (1) 

where ki and ^2 are the momentum indices corresponding to the left and right electrodes, and a is the spin index 
(conserved under tunneling); cj^^ and c^^ are creation and annihilation operators for fermions. It is then straight- 



FIG. 1: Schematic of a tunnel junction. In the left figure, the leads are the clear areas marked with 1 and 2, while the insulating 
barrier is depicted by the shaded area. A schematic of the energy levels is shown for the case of metallic electrodes. 



forward to calculate the probability per unit time (the tunneling rate) that an electron residing in the left electrode 
would tunnel from a state {ki, a} (of energy e^.j = 'h?k\/2m — epi, where epi is the Fermi level of the left electrode) 
to a state {^2, cr} (of energy = ?i^fc|/2m — eF2, where eF2 is the Fermi energy of the right electrode), by a direct 
application of the Fermi Golden rule, 

Pi^2(fci,fc2) = Y^n.M^PkjkA^- fk,), (2) 

where /fcj 2 are the Fermi functions corresponding to the states fci.2, Pk2 is the density of states (without the spin 
degrees of freedom) at and 5E is the energy difference between the states {^2, cr} and {fci, cr}, 6Ei^2 = efe2 ~ ^fci- 
The energy difference 6Ei^2 does not depend on the state fci or ^2 since it is due to electrical potentials such as an 
external battery or to the charging of the island (as we will see below in the case of single-electron transistor) which 
produce a difference in the Fermi energies of the two metals, 5Ei^2 = ePi — ^¥2-1 meaning that e/cj + epi = E/tj + ^F2, 
and thus ensuring that the total energy is conserved under tunneling. 

Replacing now the sum over states with a continuous integral, we can find the total tunneling rate from left to right 
by integrating Pi^2, weighted with the density of states pi, 

ri^2{SEi^2) = 2 J deipi{ei)Pi^2{eue2 ^ f-1 + 6Ei^2) ^ 2 J dei J de2Pi{e^)Pi^2{ei,e2)5{e2 ~ ei ^ 5Ei^2). (3) 

The last form of the equation, using the Dirac delta function (in no danger of confusion with the 5 of 5Ei^2)i reflects 
indeed the conservation of energy: electrons of energy ei + epi can get only to states 62 + eF2, if they acquire an energy 
gain of (5i?i_y2 with respect to their Fermi levels. Inelastic processes arc allowed in a more general formulation, which 
includes the degrees of freedom of the environment, in which case 5 is replaced by the so-called P- function Also, 
the factor of 2 appearing in front of the integral comes from summation over the spin index. The rate of tunneling 
from the right electrode to the left one, T2^i{5E2^i), is calculated in a similar way. The difference between these 
quantities represents the net current of charges flowing through the junction, therefore the net electrical current is 

h^2 = -e [Ti^2{5Ei^2) - T2^i{5E2^i)] - -e [ri_,2(<5-Bi->2) - T2^i{-6Ei^2)] (4) 
We then get, using the expressions Eq. (|3]) and Eq. ([2]), 

h^2 = -— I deide2|T(ei, £2)1 V(ei)p(e2)[/(ei) - /(e2)]<^(e2 - d - SE,^2), (5) 
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where we conveniently change all the momentum indices to continuous-energies functions. To make progress with Eq. 
^ , a few assumptions are needed. First of all, already embedded in the theory above is the idea that tunneling does 
not significantly change the Fermi distribution function. Nonequilibrium processes however can play an important 
role, requiring a more involved treatment using Green's functions [lo| . Next, the Fermi energy in most metals is of 
the order of a few eV^s; this is large when compared to the typical bias voltages at which such junctions are studied. 
Then, in the case of metals, the density of states can be considered constant around the corresponding Fermi energies 
p(ei) w p(0), and p{e2) ~ p(0). For superconducting electrodes this assumption does not hold: the density of states 
is divergent around the gap, and, as we will see below, this has significant effects. It is however possible (see below) 
to isolate the divergent part of the density of states as an adimensional quantity which, when multiplied with the 
normal-metal density of states around the Fermi level would give an approximate but still correct superconducting 
density of states. Finally, the tunneling matrix element T can also be taken as energy-independent, T(ei,e2) = T. 
With these assumptions, the integral over energies in Eq. ^ can be performed (see [HI, [l3 and Appendix [X]). using 
the generally useful result 



/"'fm-fi^-E)]^ I . (6) 



to give for the tunneling rates 



= -T'^'^^-^^)^(---) exp(-^WBT)-l - 

For example, at zero temperature we find ri_j.2('5i?i-^2) = (47r/7i)|rpp(eFi)p(eF2)'5£'i^2 for 6Ei^2 > and ri^2 = 
for 6Ei^2 < 0. The value of the current does not depend however on temperature. From Eq. ([7]) and Eq. Q we find 
a linear dependence in SE for the current 

Ii^2 = -—p{eFi)p{eF2)\T\^SEi^2 (8) 

If 6Ei^2 ~ ^eV, as it is in the case of a voltage-biased single junction (see Fig. [T]), we immediately recognize Ohm's 
law from Eq. ([5]), Ii^2iV) = i?j^2^y; with the quanta of resistance defined as Rk = h/e^ = 27r7i/e^ = 25.8 kJl we get 

^ = |^p(eFi)p(6F2)|Tr (9) 

-Kl2 Kk 

The theory presented above looks very reasonable until one realizes - horribile dictul - that in fact the result Eq. ^ 
predicts a rather absurd scaling of electrical resistance with the dimensions of the sample. Let us look at the densities 
of state. The three-dimensional density of states per volume (with spin states excluded) can be calculated using the 
standard text-book expression 

n^^{?^)vE. ,10) 

The Fermi energy for metals is of the order of a few eV: for example, ep^"-* = 5.53 cV, Cp"^^'' ~ 5.94 eV, ep"^'' = 11.63 
eV, e^°"^ = 7 cV. At ep = 10 eV we obtain for example g^^^^lOeV) ^ 0.75 x lO^^J-^m'^; the Fermi wavelength at 
ep = 10 eV is Ap = h/y/2rriep= 0.37 nm. 

The problem now is that in order to find /3(ep) we have to multiply the density of states per volume G^'^ with 
the volume of the electrodes. This is of course absurd, as the inverse of resistance should be linear in the junction's 
surface and not depend on the volume of the electrodes. The solution to this conundrum stems from the idea that 
the summation in Eq. ([3]) has to be done under the restriction that transversal momentum is conserved (specular 
transmission): thus the correct densities of states that appear further have to be densities of states at fixed transversal 
momentum [131 ) . Also the assumption of energy- independence of the tunneling matrix element does not necessarily 
hold true for most practical junctions. A complete theory of tunneling that takes into consideration all these special 
features has been developed over the years [ij]. However, with the trick of absorbing the density of states into the 
tunneling resistance, the theory described above is still useful, describing correctly the voltage dependence of the 
current. As we will see below, the theory can be extended to the case of superconductors. Similar tunneling theories 
can be constructed for scanning tunneling microscopy, where the densities of state refer to the tip of the instrument 
and to the surface of the sample [l^. Also, the same structure has been successfully applied to systems such as 
trapped atoms, where the tunneling can be realized by using RF fields coupled to internal hyperfine states of the 
atoms [iHl- In this case, SE is given by the detuning of the RF field with respect to the atomic frequency, and the 
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tunneling matrix element is truly constant (given by the coupling amplitude of the field with the atoms). In this 
situation, the current of particles between two hyperfine states depends indeed on the volume of the samples, as it 
should be, since the field penetrates the sample completely and tunneling between the internal states occurs not only 
at some interface but in the whole volume. This case also has its specifics however: the momentum (and not only the 
energy) is conserved. 

In the case of superconductors, the theory proceeds as above, but the density of states used has to be replaced by 
that given by the BCS theory. It turns out that this density of states consists of a factor which is the same as for 
metals, which gets multiplied by a divergent part near the gap. It is then natural to introduce a normalized density 
of states A/", which would take the value 1 for metals, 

AAnormal(e) 1, (11) 

and will be a function 

A/Bcs(e) = ^=Ji=e(|6|-A), (12) 

for BCS superconductors (A is the superconducting gap and is the Heavysidc function). 

Summarizing all the discussion above, we can write the tunneling probability, including both the metallic and 
superconducting 

oo oo 

ri^2(<5e) = Jdei j Ae2 Ni{ei)N2{e2)h{ei)[l - f2{e2)W2 - ei - 5E^^2) (13a) 

— OO — OO 

oo 

= /dei AAi(ei)AA2(ei+(5£;i_,2)/i(ei)[l-/2(ei+<5£^i^2)], (13b) 

— oo 

with A/i,2 being the normalized density of states of the electrodes, as introduced above. Furthermore, given the 
tunneling probabilities of Eq. (fT3)) . the current through a single superconducting tunnel junction without charging 
effects is given as 

h^2 = -e [ri_,2(<5-Bi^2) - T2^i{5E2^i)\ = -e [Ti^2{SE) - V2^i{-5E^^2)] ■ (14) 



B. Numerical methods 



For the discussion in this subsection we use the shorthand notation 6E = SEi^2- Now Eq. ((T^ can be solved 
analytically for some limiting cases, as we have already seen for the case of normal metal electrodes, or can be 
expressed as a combination of special functions. For example, in the case of a SIN junction, the calculation at zero 
temperature is straightforward. From Eq. (|13bp we get 

ri^2 = -^JiSEr - A?e(<5i? - Ai), (15) 

which can be seen for example in Figure [2] bottom left (the additional spike appearing there at 5E = is a singularity- 
matching peak discussed below; it shows up because the temperature is finite). In this figure, / = Ii^2 and 5E = 
SEi^2 = -eV. 

In practice however, in the general case of non-identical superconductors and finite temperature, it has to be solved 
numerically, and, due to the singularities in the BCS density of states, the numerical treatment of Eq. (jl3p is somewhat 
challenging (see Fig. [2|). We therefore want to discuss the numerical problems arising from these divergencies and 
possible ways to solve them in a bit more detail. Inserting Eq. ((T^ into Eq. (|13bp . one can easily see that the 
expression for the electron tunneling probability splits into four separate integrals, 

ri^2 = (Ti + T2 + T3 + T4) , (16) 

where the integrals T„ are given by 

min(-Ai,-A2-(5-E) 

Ti = Jde .9(e) (17a) 

— 00 
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FIG. 2: Top: Tunneling rate (left) as a function of energy difference, and I-V (right) through an asymmetric superconducting 
tunnel junction with A2 = 0.2Ai. The inset in the tunnel rate plot shows the same data with a logarithmic scale on the y-axis, 
for better emphasis of all features. Bottom: Tunneling rate (left) and I-V (right) through a symmetric superconducting tunnel 
junction. In the I-V we show the current for two different temperatures, with fesTo = 3.92 x 10~^Ai. 



T2 = -e{-Ai-A2 + SE) jdeg{e) 

A2~SE 
-A2-SE 

T3 = -e{-A,-A2-5E) J degie] 

Ai 

00 

T4 = Jde g{e) . 



(17b) 
(17c) 
(17d) 



max(Ai,A2-i5-E) 

The limits of these integrals are set by the 8-functions in the BCS density of states. The function g{e) is then 

e(e + ,5^;) 



9ie) = 



^e^-Al^{e + SE)^-A\ 



(18) 



Due to the singularities in the BCS density of states, the function g{e) is singular in each of the finite limits of the 
integrals T„. 

In most cases these singularities can be handled easily enough with the help of a simple integral transformation. 
As the function (7(e) behaves qualitatively like h(x) * x/\/x^ — c? , where h{x) is smooth on the interval of integration. 
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we can substitute x = \J ^ (? and we arc left with an integral which is numericahy easily handled: 



dx 



xh{x) 
\J x"^ — (? 



ds h{\J + c2) . 



(19) 



However, for certain values of (5i?, the singularities of the densities of states of the first and second lead will coincide, 
leading to two logarithmic divergences known as singularity matching peaks at 8E = ±(Ai + A2) and two finite steps 
at the quasiparticle threshold 5E = ±|Ai — A2I. 

Singularity-matching peaks: We will first discuss the singularity matching peaks, which are described by the in- 
tegrals Ti and T4. When 5E is close to ±(Ai + A2), the integration results can become very inaccurate, and the 
transformation done in Eq. (jl9p will not help to solve the problem. We can however separate the singularity from the 
integral by partial integration, leaving us with a finite and easily integrable rest term. We note that the general form 
of integrals Ti and T4 is 



T(c) 



dx 



h{x) 



^Jx{x + c) 



(20) 



As h{x) is proportional to the Fermi distribution, h{x) and all its derivatives exponentially tend to zero for large 
values of x. Furthermore, h{x) is analytic on the interval [0, 00). Using partial integration, we obtain 



00 

T(c) ^ \n{^^h{x = 0)-2 jdx \ll{^/^+^/^T~c)h'{x), 



(21) 



with the big advantage that the singularity itself is described by an explicit expression rather than an integral. The 
integral in the second term is finite and easily computed. The major nuisance here is to calculate the derivative of 
the function hix). 

For completeness, it is also possible to derive a complete expansion of T(c) with the help of repeated partial 
integration. 



T(c) = lnQ)Mx = 0)-f:igj^c"ln(c)/.(")(. = 0), 

^ ^ 71—1 ^ ' 



(22) 



Here we denote by h^'^\x) the nth derivative of h{x) with respect to x (for a derivation, see Appendix IB|) . 
Using the recipe of Eq. (j21[) on the integrals Ti and T4, we obtain 



Ti = Si - 2 J deln 
max(Ai,A2+<5-E) 



\/e- Ai + y^e-SE - A2 



de 



e(e - ,5^;) 



v/(e + Ai)(e-(5S + A2) 



[l-fi{e)]h{e-6E) 



(23a) 



T4 = '5'4 - 2y deln 
max(Ai,A2-i5B) 



y/e~ Ai + ^/e + SE - A2 
e(e + SE) 



de 



V(e + Ai)(e + (5S + A2) 



fi{e)[l- f2{e + 5E)] 



(23b) 



where the terms Si and ^4 describe the singular points, with different expressions if the singularity is approached 
from above or below, 



SiidE > Ai - A2) 
Si{SE < Ai - A2) 



A2(A2 + SE)h{-A2 - SE)f2{A2) 
^2A2(Ai+A2 + 5i;) 

AiiAi-5E)hi-Ai)f2{Ai-6E) 
^2Ai{Ai+A2~SE) 



ln(-Ai + A2 + 6E) 
ln(Ai - A2 - 6E) 



(24a) 
(24b) 
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SiiSE > A2 - Ai) = - 
S^iSE < A2 - Ai) = - 



Ai(Ai+(5£;)/i(Ai)/2(-Ai-^^) 

A2(A2 - 6E)h{A2 - SE)f2{-A2) 
^2A2{A2+Ai-6E) 



ln{-A2 + Ai+ SE) 
ln(A2 - Ai - 6E) . 



(24c) 
(24d) 



In the case of a symmetric junction, where Ai = A2 = A, only one singularity matching peak appears in the tunneling 
rates (jl3p . In this case the singularity matching peak disappears in the current through the junction. With the derived 
expressions we can easily understand this, as the expressions 5*1 and ^4 cancel out when inserted into the expression 
for the current (jl4p . In this special case the current can therefore be computed using only the respective second terms 
in Eqs. 

Quasiparticle threshold: The quasiparticle threshold is described by the integrals T2 and T3 in Eqs. PT]) . We can 
find an explicit expression for the height of these jumps with the help of a partial integration and taking the limit 
6E — >• ±(Ai + A2). Both T2 and T3 are of the general form 



c 

T(c) = Jdx 



h{x) 



(25) 



where both h(x) and h! (x) are analytic functions on the interval 
limit c — >■ by partial integration. 



Then we can find the solution of T(c) in the 



Tc — -^\h(c) ^ h[^—c)\-- /dxarctan 



7r/i(0) + , for c ^ 



h!{x) 



Using this method, we find that the height of the step at the quasiparticle threshold voltage is 



2e2i?i2 



V'AiA2/i(-Ai)/2(-A2) 



for bE = - {Ax ^ A2), and 



VAiA2/i(+Ai)/2(+A2) 



(26) 



(27a) 



(27b) 



for5£;=:+(Ai+A2). 

We note here that at zero temperature the height of the first step is {-k l2e?R\2)\J A\A2 and the second one is zero 
(due to the fact that Fermi functions at positive energies are zero). The first peak is always more prominent and 
positive, while the second threshold is observable only at finite temperatures and it appears as a negative step (see 
the inset of Figure . 

Dynes parametrization: As one can see from the above discussion, it is rather cumbersome to treat superconducting 
tunnel junctions in the BCS picture. The expressions for the tunneling probabilities of Eq. p3p will also appear in the 
treatment of superconducting SETs, where the strong features discussed in this chapter will cause further numerical 
difficulties, especially when computing conductances. It is therefore convenient to somewhat smoothen the jumps and 
singularities of the tunneling probabilities. 

One method to achieve this is to replace the (^-distribution in Eq. (|13a[) by a Gaussian function of finite width. 
(This is not only a mathematical trick but it corresponds, in the so-called P{E) theory [ll| to the case of a high 
impedance environment at finite temperature). One is then left to solve a two-dimensional integral with non-trivial 
integration limits. 

A more convenient and at the same time physically more justifiable solution to the problem is to introduce a 
broadening parameter into the density of states (fT2)) . 



Dynes 



E 



IT] 



^{E + IT^Y + A2 



(28) 



The parameter rj is called Dynes parameter [17[ and it accounts for the finite life-time of quasiparticles in the super- 
conductor. The symbol 5R stands for the real part of the expression in curly brackets. This broadening of the density 
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of states leads to nonzero currents flowing through the tunnel junction at voltages smaller than A/e (say for a NIS 
junction), which are often seen in experiments. These sub-gap currents are believed to be due to additional states 
formed inside the gap, allowing for the opening of new conduction channels. The origin of these currents is still not 
clarified in the literature. One should realize that these subgap currents are distinct from the currents due solely to 
the existence of a finite temperature (which causes subgap excitations) in a superconductor. Such temperature effects 
are already included in the formulation presented above. 

It has been demonstrated [l^ that these additional currents due to states formed inside the gap are Coulomb- 
blockaded as well (as it could be expected physically). Moreover, they can be even used for Coulomb-blockade 
thermometry. It should be noted also that an intrinsic dependence of the gap on temperature and magnetic field exists, 
as given by standard BCS theory. By analyzing the change in the zero-bias conductance when either temperature or 
magnetic field is increased, it is possible to determine the critical temperature and the critical magnetic field of the 
island [l^ . A detailed description of another fitting procedure which uses the Dynes form for the density of states is 
also given in e.g. [IJ]. Through replacing the density of states in Eq. ([T^ with Eq. (^5)) the tunneling probabilities 
in Eq. ([T^ become finite and continuous for all values of 5E and integration takes place on the whole real axis. The 
integral (jl3b[) is still not easy to solve with high precision, but compared to the procedure that follows when using 
Eq. (|12p. the numerical treatment becomes certainly more simple. In the following sections we will exclusively use 
Eq. (pS)) as density of states, with a small but finite value for 1]. 

III. TWO-JUNCTION DEVICES 

Two-junction arrays are very useful devices for a variety of applications. They are denoted usually as a het- 
erostructure - for example, SINIS, NININ, etc., showing the phase of the electrons (metallic or superconducting) at 
the operating temperature. If a gate is capacitively coupled with the middle electrode (called island) by a coupling 
capacitance Cg, then the device is usually referred to as a single electron transistor (SET). If some or all of the leads 
in the SET are superconducting, the device will have characteristics that stem altogether from tunneling, charging 
effects, and superconducting gap. The fabrication of such devices is possible due to modern nanolithography tech- 
niques, as mentioned in the previous section, using which one can easily fabricate small junctions and small metallic 
or superconducting grains, resulting in non-negligible charging energies. For example, the area S of the junctions 
obtained by shadow evaporation techniques can be of the order of 100 nm^ or even less. The insulating layer is usually 
of the order of d =1 nm (usually 0.5-5 nm) in thickness, and the dielectric constant of the oxide is of the order of 10. 
Using the formula for the capacitance C = eoerS/d, we can readily obtain capacitances of the order of C « 10~^^F 
(fcmtofarads). The energy scale corresponding to this capacitance is e^/2C sa 100 fiV (or a temperature of the order 
of IK). Thus, charging effects can become observable at temperatures that can be achieved by standard cryogenic 
devices, such as dilution refrigerators. Effort has been put in recent times in reducing this capacitance even further, 
so that charging effects would be visible at even higher temperatures (preferably at room temperature [1^). This is 
where novel fabrication methods, based on advances in nanoscience, are essential. 

A. Master equation 

The Coulomb blockade phenomena discussed here are essentially classical. The tunnel resistance is assumed to 
have a value much larger than the quanta of resistance Rk = h/e^ . This makes the number of electrons on an island 
a well-defined classical quantity, albeit a discret one. The Coulomb blockade effects then stem from the fact that 
the biasing potential can be varied continuosly while the potential given by the product between the charge and the 
capacitance can change only in steps. This gives rise to steps in the current and Coulomb blockade oscillations in 
the conductance. Except for tunneling, this phenomenon is indeed a classical one: the energy level separation of a 
small metallic island is still much smaller than the thermal energy ksT, therefore the spectrum can be treated as 
continuous. This should be contrasted to the case of tunneling in semiconductor nanostructures, where the spectrum 
is discrete (see [2l|). 

In the case of single electron transistors, the relevant charging energy is given by the effective capacitance of the 
island to the ground: this is denoted in the following, using an already standard notation, by Cs, expressed as 
Ce = Cl + Cr + Cg. Here Cl is the capacitance of the left junction, Cr is the capacitance of the right junction, 
and Cg is the capacitance between the island and the gate electrode. The associated charging energy of the island, 
corresponding to one electron of extra charge, is 



(29) 
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FIG. 3: Schematic of an SET. The top picture 
shows the left (L) and right (R) leads in contact 
with the island (I) via a tunneling junction, each. 
The bottom picture shows the electrical equiv- 
alent circuit, where Vl and Vb. are the voltages 
applied to the left and right lead, respectively, 
while Vg is the gate voltage, that effectively al- 
ters the island charge. 



The gate capacitance Cg is usually neglected or it can formally be distributed between the junction capacitances, 
leaving Ce « Cl + Cr. The relevant electro-static energy of the island with n excess electrons and gated with the 
potential Vg turns out to be 

1 2 
Echin,Qg) = 26^^'^^^^^-' 

EE Ec{n + Ugf , 

where Ug = qg/e and Qg = CgVg + CrVr + ClVl- We want now to calculate the energy changes SEi^^i{n), SEi^i^{n), 
(5i?R_s.i(ri), and 5Ei^B. {n), (including the work done by the sources) associated with electrons tunneling onto and 
off the island). The general idea is that, when such processes occur, n is increased (or decreased), respectively, by 
one, and the energy difference between the electron states before and after the tunneling event will be given by 
Ech{n -f 1) — Ech{n) (or Ech{n) — Ech{n — 1)) plus the work done by the voltage source. 

We give now a detailed derivation of these results. Kirchhoff's laws for the circuit presented in Eq.(|3]) can be written 
as: 

)g, (30a) 

(30b) 

(30c) 
(30d) 

where Ql, Qk a-nd Qg are the charges corresponding to capacitors Cl, Cr and Cg respectively. Combining these 
equations, we find 

= ~^[C^yK + Ciyi^-Vg{Ci^ + C^)+ne] (31a) 
Qi.{n) = ~^[C^yK + CgVg-Vi^{Cg + C^)+ne] (31b) 
QrW = -^[Ci.Vi^ + CgVg~VK{Cg + Ci^)+ne] (31c) 
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The total electrostatic energy of the system stored in the capacitors is given by 

Using the expressions Eqs. ([31]) one can verify immediately an interesting property of E^'^^\n), namely that it does 
not contain any terms linear in n, 

E^°^\n) — — — [-terms depending only on voltages. (33) 

Consider now a process by which one electron tunnels from the left electrode to the island, accompanied by a 
re-arrangement of the charges so that electrostatic equilibrium is reached corresponding to a state with n + 1 extra 
electrons on the island. The work done by the sources is then given by 

W^L-yi = [QLin + 1) - QL{n) - e] + Vk [Qnin + 1) - QkH] + Vg [Qg{n + 1) - Qg{n)] . (34) 

Let us look a bit in slow motion at what is in fact happening here: as the electron tunnels through the left island, 
to ensure the neutrality of the conductor connecting the source Vl to the left island, another electron must be pulled 
through the source, which requires the energy — cVl. As the electron arrives on the island, the system is no longer 
in the electrostatic equilibrium ensured by Eqs. ([50]) and, to reach the new equilibrium with n-\- 1 electrons on the 
island, the charges Q-L{n -I- 1) — (5l("-), Qr('t- + 1) ~ Qr("'), and Qg{n -|- 1) — Qg{n) are being transferred through the 
corresponding sources. The total energy change during this process, which includes the work done by the sources, is 

SEi^^i = Wi^^i + E^'''\n) ~ E'-'^'^n + 1). (35) 

In other words, when an electron tunnels from a state fcL of the left electrode to a state qi of the island, the energy 
is conserved (due to the Dirac delta appearing in the transition rate) in the following way: eqj ~ e^^ + SEj^^i, or 
eg,+E'-''^\n + l) =€k^+WL^i+E<^''^^n). There are altogether four different possible tunneling events, two increasing 
and two decreasing the island charge: 

SEi^^iin) = -cVl - 2Ec {n + Ug + 1/2) , (36a) 

SEi^i^in) = +eVi + 2Ec {n + Ug- 1/2) , (36b) 

SEi^nin) = +eVR + 2Ec {n + Ug - 1/2) , (36c) 

SEK-^i{n) = -cVb. -2Ec{n + ng + 1/2) . (36d) 

Another handy way to get these changes in energy is to introduce the free energy of the system as the Legendre 
transform of the electrostatic energy Eq. (j32[) 



E,i,{n) = E^°'\n) - gL(n)VL - QR{n)VR - Qg{n)Vg. (37) 
Then, using Eqs. ((3T|) we get (up to a constant) 

EM = T^ine + qg)\ (38) 

and the expressions Eqs. p6p can be put in a physically more transparent form 

6EL^i{n) = -cVl + E,i,{n) ^ E,h{n + I) , (39a) 

6Ei^i^{n) = +eVL + E,i,{n) - E,i,{n - I) , (39b) 

SEi^uin) = +eVR + E,,,{n) - E,h{n - I) , (39c) 

SEn-^i{n) = -cVr + E,hin) - E,h{n + I) . (39d) 

We see that Ech{n) has the meaning of an effective charging energy of the island, including the work done by the 
sources. 

The probability to tunnel onto or off the island through either the left or right junction can be obtained by inserting 
the respective energy difference of Eq. p6p into Eq. (fT5)) . We define the probabilities for the four different tunneling 
events as function of the excess charge n on the island as 

rL^i(n) = T[5E^^i{n)] (40a) 

ri^L(n) = T[5Ei^i,{n)] (40b) 

ri^R(n) = T[5Ei^nin)] (40c) 

rR^i(n) = T[5En^i{n)]. (40d) 
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To continue our analysis, we now notice that, due to the fact that the number of electrons on the island is a well- 
defined number, we can talk about classical states of the system and index them with the excess number of electrons 
on the island, n. There will then be transitions between these states, namely, with the notations in Eqs. (|40)) . the rates 
corresponding to the the island's charge number n being raised or lowered by one (T^ and respectively) become. 



rL^i(n) - 



rR,^i(n) . 



(41a) 
(41b) 



We now want to calculate the rate at which the probability p{n; t) that the system is in state n changes. There are 
three contributions to this rate: 

• If the system is already in the state n, it can make transitions to the state n + 1 or to the state n — 1, with 
transition rates r-|-(n) and r^(n); 

• The system can get to the state n from state n — 1 with transition rate r-|-(n — 1); 

• The system can get to the state n from state n + 1 with transition rate r^(n + 1). 
This results in a master equation for p{n] t), 

dp{n; t) 



dt 



-[T^{n) + T^{n)]p{n- 1) +T^{n+ l)p{n + 1; t) + T^{n - l)p{n - 1; t). 



(42) 



This equation describes a classical Markovian process: indeed, embedded in this treatment is the idea that electrons 
do not have any memory while tunneling. The history of the electrons being at previous times in other states is erased 
and therefore there is no time-dependence in the transition rates Eqs. ([1T|) . which depend only on the "present" state 
n. 

We now want to characterize the stationary state, defined by the condition dp{n;t)/dt ~ 0. From Eq. ([^^ we can 
verify immediately that the stationary state is given by a time-independent p{n) satisfying 



p{n)r^{n)=p{n + l)r^{n + l), 

which is equivalent to 

p(n)r^(n) = p{n - l)T^{n - 1), 
Using now the normalization condition for probabilities, 

oo 



(43a) 
(43b) 

(44) 



wc find 



p{n) 



U rt(0 

i — — oo 



n r,(*) 

_i—n-\-l 



m— — oo 



?n — 1 

n r^(0 

i——oo 



n r,(z) 

i—m-\-l 



(45) 



From Eqs. fj43p another useful property results immediately, 

oo oo 



(46) 



n— — oo 



n— — oo 



Using now Eq. P5|) and the definitions in Eqs. (jU]) it is easy to check that there is no charge accumulation on the 
island, that is, the currents fiowing through the left junction and the right junction are equal. 



Y P^") [rL^i("-) - Ti^L{n)] = -e E p{n) \ri^R{n) - rR^i(n)] 



(47) 



n— — oo 



n— — oo 



This formula is used in the following to calculate the currents and the conductivities. 

We point out here that the simplifications leading to Eqs. (|44ll47p are valid only under the assumption of steady- 
state conditions. In cases in which biasing is done by RF fields - such as turnstile operation or RF cooling - this 
assumption could break and the treatment above is insufficient. 
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FIG. 4: The current (left) and the conductance (right) of a Coulomb blockade thermometer at different temperatures. For 
temperatures T < Ec one can clearly see the steps in the current and the oscillations in the conductance. For T > Ec these 
steps/osillations smear out and the Coulomb dip forms. The inset in the conductance plot shows the Coulomb dip curves 
(T > Ec) in a voltage range of [— 100-Ec, 100-Ec] to make the typical shape of the dip also visible for higher temperatures. 



B. Coulomb blockade thermometry 

A very useful application of double junction structures with normal-metal electrodes (NININ) is for temperature 
measurements. Consider such a structure, with charging energy of the island denoted as before by Eq = e^/2Cs and 
at temperature (assuming equilibrium) T. Clear Coulomb-blockade effects can be seen if the temperature is much 
smaller than the charging energy: the current is zero for voltages lower then e/Cs, as it takes that much energy to 
add an electron on the island. Plateaus are then formed at values of the bias voltage corresponding to multiples of 
e/C (see Figure U] left). The conductance curves then show so-called Coulomb-blockade oscillations (see Figure U]) 
right). Gate voltages can be also included in this analysis, in which case it can be shown that, in the gate voltage ~ 
bias voltage space the stable states (corresponding to different numbers of electrons on the island) form a diamond 
pattern |23l |. 

We now examine the effect of temperature: intuitively, as the temperature is increased, one expects that the features 
visible at low temperatures will be smeared out and eventually washed away. Surprisingly however, if the temperature 
is larger than the charging energy, one clear feature survives in the conductance, namely a dip at low bias voltages. 
It can be shown within the orthodox theory presented above, that, in the limit fcsT ^ Ec, the full width of 
this Coulomb blockade dip (measured half-way between the dip in the conductance at zero bias and the plateau at 
V — > ±oo) is given by 

Tjr^ « 5.439, (48) 

and the relative change of conductance at zero bias voltage is 

G{V ±oo) - G{V = 0) Ec 



G{V ±oo) 3KbT 



(49) 



It is straigthforward to show that G{V — ±oo) = (2i?)~^, meaning that at large bias voltage the charging energy of 
the island docs not have any effect and the conductance of the NININ array becomes that of two junctions of resistance 
R in scries. Eq. (|^5)) enables the use of such structures as primary thermometers: from a simple measurement of 
conductance it is possible to extract V1/2 without any need for calibration. Eq. (j49p can be used to determine Ec- 

What happens now if the electrodes become superconducting? In Figure [5] left we present the behavior of the 
Coulomb blockade dip for a SINIS structure as the gap of the electrodes is increased. One can see how the supercon- 
ducting gap is introducing additional features into the conductance curve which, as the gap is gradually increased, 
finally distort its shape so much that the CBT dip is very soon not distinguishable anymore. The origin of these 
features is described in detail in Section IIIIDI 
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FIG. 5: Left: The conductance dip in a superconducting CBT at a temperature of T = 5Ec, with a finite gap A/ = Al ~ 0-lEc 
(red color). For comparison we overlayed the grapli witli the corresponding curve for the all- metallic case (black line). The CB 
dip is still clearly visible, but the features due to superconductivity are clearly imprinted on the curve, according to Eq. (|51|) . 
Right: The conductance dip in a superconducting CBT, with a growing superconducting gap from to O.SEc- Here we can 
see how the dip is gradually smeared out by the features that arise due to the superconducting nature of the island and leads. 



C. SINIS 



A supcrconductor-insulator-normal metal-insulator-supcrconducting (SINIS) device can be seen as an SET without 
gate, with superconducting leads and a normal metal island. SINIS are very useful devices, used as local coolers 0], 
secondary thermometers [24 | , and electron pumps [25| . Usually the charging energy in a SINIS is so small that it can 
be neglected. However, in some applications, for instance as thermometers, the SINIS dimensions may become small 
enough to make charging effects become observable. 

Figure [ni shows the I-V-cmves and conductances of a SINIS structure at different operating temperatures for both 
zero and finite charging energy. Without charging energy, the SINIS reproduces the I-V curve of a single NIS junction, 
with all features being at twice the voltage, as we expect due to the fact that only half of the voltage applied to the 
SINIS drops over a single junction. With a finite charging energy, the quasiparticle threshold is pushed up from 
V = 2A/e to (2A + 2Ec)/e, and is repeated at AEc intervals with decreasing amplitude. This is visible especially 
well in the conductance plots in Fig. [6] At operating temperatures that are comparable to or larger than the charging 
energy, all features are smeared out. In Fig. [5] (for the case that A = 220/i eV and thus Ec = 33/^ eV), the curve 
with still barely visible oscillations conductance would correspond to a temperature of 2kBTo « 17/i eV, while in the 
curve at the next higher temperature of bksTo « 43/x eV the oscillations are smeared out entirely. 



D. Superconducting SETs 

When also the island becomes superconducting, we differentiate between two cases, the symmetric SSET, where 
the gaps (denoted by A) of both leads and island are identical, and the asymmetric SSET, where the gap in the island 
differs from the gaps in the leads (denoted by Aj and Al, respectively). For simplicity of the further discussion, we 
introduce two new notations, namely the sum of the lead and island gaps. 

As = Al + Ai , (50a) 

and the difference between the gaps, 

A5 = |Al-Ai|. (50b) 

Figure [7] shows the I-V characteristics of an asymmetric SET with Aj — 0.2Al and Ec — for different tempera- 
tures. Compared to the I-V of a SINIS structure, we can see the appearance of a new feature, the singularity-matching 
peak, which apears at voltages = ±2As/e. As mentioned above, these peaks appear because of the overlap of 

two infinite density of states appearing under the integral giving the tunneling rates. Singularity- matching peaks in 
superconducting single-electron transistors have been investigated theoretically previously in [26[ and first observed 
experimentally by (231 and [g^- In NbAlNb structures subgap features at 2(ANb — Aai) have been observed in [29} . 
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FIG. 6: Top: Current (left) and conductance (right) of a SINIS structure as function of applied voltage at different operating 
temperatures. The charging energy of the system is set to zero. Bottom: Current (left) and conductance (right) through a 
SINIS structure with finite charging energy {Ec = 0.15 A). The offset of the quasiparticle threshold is 2Ec and the periodicity 
of its repetition is AEc- In all plots the temperature is ksTo = 3.92 x 10~^A. 



Also in this figure we see the usual feature due to quasiparticle threshold at y — ±2As/e. Finally, a third peak 
appears at V — 2Ai^/e, where the upper and lower singularity in the density of states of the left and right lead, 
respectively, arc aligned. The right plot in Figure [7] shows the conductance for the same values. The conductance 
plot is perhaps not as descriptive as the one of a SINIS, but it still serves two purposes: 1) it emphasizes very well 
all sharp features and thus makes it easy to determine their position and, 2) it allows to easily distinguish between 
singularity matching peak- and quasiparticle threshold-like features, as the latter only produces an upward spike, 
while the former produces a double spike, with the second spike going downwards. This is especially useful when 
discussing SSETs with charging energy, where one can find a large amount of weak features. 

At finite charging energies, we see a similar repetition of features as one sees in SINIS I-Vs with Ec > 0. These 
features are a combination of the usual coulomb blockade features in normal metal SETs and a series of peaks and 
steps. The origin of the peaks are the singularity matching peaks in the tunneling rates (jl3p . Their positions are 
easily calculated by equating the energy differences of Eq. (j36|) with the positions of the singularity matching peaks 
ibA^, and solving for the applied voltage. 

To get a feeling of how these features enter the I-V-cmve, we first look at the special case of two identical junctions, 
where Cl = Cr = Cs/2. In this setup, we have that Vl = —Vr ~ V/2. For further simplification, we also set the 
gate charge equal to zero, i.e. Ug = 0. If we now for instance look at the series of peaks that correspond to SE^^j, 
we find ±As = ey(™V2 - 2Ec{n + 1/2), and therefore eF^"™) = ±2A5 + 2Ec{2n + 1). Altogether we find that the 
bias voltages at which peaks occur (at Ug ~ 0) are 

yi^-)=±^±^(2n + l). (51) 
e e 

Similarly, the I~V shows a series of positive and negative steps. These steps originate from the quasiparticle threshold 
in the tunneling rates (|13p . We can find them in a similar manner as the peaks, by equating Eqs. p6p with itAs and 
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FIG. 7: Current (left) and conductance (right) of an asymmetric SET as function of applied voltage at different operating 
temperatures. The charging energy of the system is set to zero, fcsTo — 3.92 x 10^'^ Al, and Ai = 0.15Al, where the subscripts 
I and L denote island and lead, respectively. 
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FIG. 8: Current (left) and conductance (right) of a symmetric SET as function of applied voltage at different operating 
temperatures. The gaps are chosen such that Ai — Al = A, and the charging energy of the system is set to Ec = 0.15A. As 
before, fcsTo = 3.92 x 10~'^A. In both plots the top most scale denotes offsets from the value of eV — 4A. 



find 

y^(qp)=±^±^(2n + l). (52) 

e e 

Equation (j52p corresponds to positive steps for and to negative steps for —As. The negative steps however, 

vifhich are already small in the tunneling rates, are difficult to observe in practice. 

Figures [HI m and [TU] shovif the I~V and conductance of SSETs vifith a finite charging energy of Ec — O.lSAi, with 
ratios between island and lead gaps of Aj/Al = 1, 0.5 and 0.2, respectively. In all figures one can see the distribution 
of peaks and steps as described in Eqs. (jSip and (1521) . with growing complexity. In the I-V curves, especially in Figure 
IHl the features are difficult to distinguish. Therefore the conductance plots become a powerful tool in identifying the 
position of peaks and steps, and in distinguishing between them (single and double spikes in conductance correspond 
to steps and peaks, respectively, in current). 



E. Effect of the gate 



To study how the gate charge affects the IV-curves of an SSET, it is best to do some simplifications again, to 
somewhat limit the amount of confusing details. In Fig. [TT] we have reproduced the first figure of where the IV 
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FIG. 9: Current (left) and conductance (right) of an asymmetric SET as function of applied voltage at different operating 
temperatures. The gaps are chosen such that Ai = 0.5Al, and the charging energy of the system is set to Ec ~ 0.15 A. As 
before, fcsTo = 3.92 x 10~^Al. In both plots the top most scale denotes offsets from the value of eV = 2Ae, while the bottom 
most scale denotes offsets from the value of eV = 2 As- The three unmarked peaks correspond to —2 As + lOEc, 14i?c and 
18Ec, respectively. 




FIG. 10: Current (left) and conductance (right) of an asymmetric SET as function of applied voltage at different operating 
temperatures. The gaps are chosen such that Ai — 0.2Al, the charging energy of the system is set to Ec = 0.15A, and 
ksTo — 3.92 X 10~'^Al. In both plots the top most scale denotes offsets from the value of eV = 2Ae, while the bottom most 
scale denotes offsets from the value of eV = 2(5A.The unmarke peaks correspond to —2 As + lAEc- 



of a symmetric SSET was studied. The figure shows a set of IV curves through the SSET as function of bias voltage 
for different values of the gate charge. As one can see, also here appears a series of peaks and steps, but their position 
is shifted as the gate charge is altered. Again the position of the features is easily calculated by equating the energy 
differences (for the most general case, use ((39)) ) to itA^ for peaks, and ±As for steps. As in this case A^- = 0, the 
equation for the position of the peaks is somewhat simplified, as here the derived voltages will not depend on the 
superconducting gap. As in [l^l, the positions of the peaks in the system at hand (where possibly Cl ^ Cr) are given 
by 



{L,R} 



1 , 



(53) 



Here we defined n'., 



Qa/e, with eng the gate charge and Qq the background charge of the island. Accordingly, 



the positions of the positive and negative steps are given by 



T/qp 

{L,R},n 



{L,R} 



ClCr 



, 2ACs 1 . 



(54) 
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FIG. 11: The effect of the gate charge. Here A — 2.6Ec and T = O.SEc- The curves are offset by the respective gate charge 
for better distinguishabihty. 



Here the positive steps are obtained for the + sign in from of the fraction 2ACs/e^, and the negative steps are 
obtained for a — sign in front of the same fraction. The negative steps can be observed only at higher temperatures 
- see the comments on the height of the quasiparticle threshold, Eq. ([?7)) and (P7|) : only the positive steps are visible 
in Figure (fTTj) . 

It is easy to verify that Eqs. ([51]) and ([5^ are another case of the expressions derived above Eqs. (|54ll53p : these sets 
of equations become identical when using the same assumptions (identical junctions and 71^=0), and remembering 
that here = and A/ = Al = = A. 



F. Cooling 

Tunneling of electrons through junctions results in energy being transferred between the electrodes, as the higher 
energetic electrons are extracted or pumped from or into one of the electrodes. In the case of two-junction structures 
such as SETs, SINIS, etc., this results in the heating and cooling of the middle electrode (the island). This is due to 
the fact that, although the currents through the two junctions are the same, as it should be due to the conservation 
of charge (the charge of the island has to remain constant) the energy transfer (heat) does not need to satisfy such a 
conservation law. 

It is straightforward to calculate the energy transferred per unit time out of the island (cooling power) as a function 
of the applied voltage, 

oo 

= J2 P^"^^ [Pi^Lin) - PL^iin) + Pi^R(n) - PR^i(n)] . (55) 

n— — oo 

Here the powers Pij^i{n), Pi^i^{n), PR_yi(n), and Pi^Y{ {n) correspond to energy being transferred onto/off the island 
by tunneling into the left and right electrodes, at a given number of excess electrons n on the island. They can be 
calculated by multiplying the energy carried by each tunneling electron to the probability of tunneling per unit time 
(given by the density of states and Fermi factors), and then summing over energy states (see also Eqs. (|13ap . (jl3bp ). 

oo 

Pi^xin) - [dE ENi{E)Nx[E^5Ei^x{n)]fi{E)[l- fx{E^SEi^x{n))] (56a) 

e^Rix J 

— oo 
oo 

Px->i{n) = fdE[E + 6Ex^i{n)]Ni[E + SEx^iin)]NxiE)[l-fi{E + 6Ex^iin))]fxiE), (56b) 

e^itx/ J 



where X = L, R. 

In Figure [T^ we present a comparison between the cooling powers on the island Eg. ([55]) for SINIS (upper plots) and 
superconducting SET structures (lower plots). We observe that cooling exists only in an interval of bias voltages not 
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FIG. 12: Cooling powers in different double-island devices. Top left: SINIS structure with Ec = at different temperatures. 
Top right: SINIS at fixed temperature STo with different charging energies. Bottom left: SET with 5To and Ec = for different 
gap sizes Ai < Al of the island. Bottom right: SET with Ai = 0.2Al and Ec = 0.15Al at different temperatures. In all plots 
ksTo = 3.92 X 10"^Al. 



much larger than the quasiparticle threshold at 2A/e (we take identical superconducting leads with gap A). As seen 
from the left upper plot, as the temperature is lowered the process becomes much less efficient. The general effect of 
the charging energy is that it is detrimental to cooling (right plot, upper part). In the case of superconducting SET's 
(lower plots, Figure [T2|) the existence of singularity-matching peaks produces relatively sharp spikes in the cooling 
power. Essentially, the island is cooled down by BCS quasiparticles [301 ■ The left figure shows these features for the 
case of negligible charging energy (a SISIS structure with a relatively large middle electrode having a different gap 
than the leads). An interesting feature is also that the range of voltages over which cooling occurs is extended, due 
to the fact that now the quasiparticle threshold is at values 2Al + 2Ai. Finally, cooling in SET structures with finite 
charging energy is shown in the lower-right plot of Figure 1121 



IV. CONCLUSIONS 



We have presented the theory of tunneling in metallic and superconducting two-junction arrays such as single 
electron transistors, together with a number of applications. All three energy scales, the charging energy, the super- 
conducting gap, and the temperature, are considered and their role is thoroughly discussed. For example, we examined 
how a finite superconducting gap affects the Coulomb-blockade based thermometry, the effect of singularity-matching 
peaks in the current-voltage and conductance-voltage characteristics of superconducting single-electron transistors, 
and we looked at the effect of charging energy in cooling devices. With the development of the field of nanotcchnol- 
ogy, such devices could emerge as very useful tools for high-precision measurements of nano-structured materials and 
objects at low temperatures. 
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Appendix A: A useful integral 

For completeness, we give a full derivation of the result Eq. ^ (see e.g. [HI). The Fermi functions are given by 

We now introduce the function 

oo 

g{x)^ J de[fie)-f{e + x)] (A2) 

— OO 

Obviously ^(O) — 0. The first derivative of g gives 

oo 

- d(e + x)— — ■ — r- = -/(oo) + /(-oo) = 1. (A3) 



dx J d{e + x) 

— oo 

This means that g{x) = x, so we have the result 

de[f{e)-f{e + x)]^x. (A4) 



Now. integrals over expressions of the type f{E)[l — f{E + x)], which appear in the theory of tunneling due to Pauli 
exclusion principle, can be solved by noticing that 

where /3 = l/ksT. Now, using Eq. (jA4p we get immediately 



Appendix B: Derivation of some integral expansions used 



Generic procedure: In the following we want to derive the expansion of Eq. (j22p . To get a feeling for the procedure, 
we will first consider a somewhat easier and more common integral, namely 

oo 

T(c) = [dx -^h{x), (Bl) 
J x + c 



which we assume to not be solvable analytically. For the function h{x) we assume that it is well behaved for all 
x e [0, oo) and that h{x) behaves as a exp(— x) for large x. However, T(c) is logarithmically divergent for c = 0. We 
can separate T(c) into a divergent and an non-divergent part by performing a partial integration, 

oo 

T(c) = \n{l/c)h{0)~ Jdx \n{x + c)h'{x). (B2) 



Here the first term is an explicit expression that describes the divergence around c = and the second term is a finite 
integral that needs to be computed numerically (we use the notation h'{x) = dh{x)/dx). Unfortunately, the function 
over the integral in the second term of Eq. (|B2|) has no Taylor expansion in a; + c — >■ 0, making any further analysis 
of the integral difficult. It is however possible to find an expansion by doing repeated partial integration. If we define 
the function Fn{x) such that 

-^K(x) = ln(x), (B3) 
dx^ 
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it is easy to show that 



Fnix) = i-x"[ln(x)-i/„], (B4) 



where = X]fc=i ^/^ harmonic number (we could actuaUy still add an arbitrary polynomial of nth order to 

Fn, but due to the repeated partial integrations all its coefficients cancel out and we therefore set them equal zero). 
These functions have the properties Ff){x) = ln(a;), F„(x = 0) = 0, (V) n > and F„-yoo{x) — > 0, (V) a; > 0. With 
this we can write T(c) as 

OO _j 

T(c) = ln(l/c)/i(0)-^(-l)"-c"[ln(c)-i/„]/i(")(0), (B5) 

where by h^"-\x) we denote the nth derivative of h{x), d"'h{x) / dx"" . 

Application: In Section |TT1 we have encountered slightly more complicated integrals, which were of the form 



T(c) = Ax —===h{x), (B6) 
J \J x{x + cj 

with a logarithmic singularity in c = 0. For this reason we would like to expand T'(c) after the recipe that was 
used when expanding the previous integral, Eq. (jBl|) . The only complication we encounter is that we need to find a 



function G„ {x) with dG^ (x) /dx ~ 1/ \/x{x + a) and 

— G„(a:) = G,{x). (B7) 

It is easy to verify that the function Go is given by 

Go(x) ^ 21n(^/^+^/^T^). (B8) 
The general function G„ can be found with the trial solution 



Gn{x) = Pn{x) ln(V^ +Vx + c) + Q„ix)^/x{x + c) , (B9) 

where Pn{x) and Qn{x) arc polynomials. By inserting Eq. (jB9l) into Eq. (|B7p . we find that P„ and Qn must fulfill 
the relations 

= (BIO) 

and 

Fn{x) + i2x + c)Qnix)+2xix + c)^^^ = 2xix + c)Qn-iix) . (Bll) 
These relations can finally shown to be fulfilled if 

n 

Pn{x) = ^^an-fcx'^ (B12a) 

fc=0 

Tl-1 

Q„(x) = (B12b) 



where 



fc=0 
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Although these polynomials look somewhat complicated, only the constant term in P„ (x) will enter the expansion of 
T(c), 



T(c) = -^G„(0)/i(")(0) 

oo 

= -^a„ln(^^)/l(")(0) 

Tl = 

= ln(^)MO)-|:^c"ln(c)/.(")(0), 



(B14) 



where the first term of Eq. (jB14|) contains the convergence for c — > 0. 

The convergence of expression (|B14p is somewhat tricky to show as it depends on the values of the function h{x) 
and its derivatives at a: = 0. With the ratio test of convergence for infinite series we find, writing Eq. (|B14[) as J^n '^n, 
that 



dn+l 



(2n+l)(2n + 2) h^'^+^\0) 
-c- 



4(n+ 1)3 

-^ln(/i(")(a;)) 
dx 



/i(")(0) 



for n — ?> oo , 



(B15) 



x=0 



which is smaller than one if (rf/da;)ln(/i(")(x))| x=a increases slower than n. For example, if h{x) = cxp(— ax), we 
would have |a„+i/a„| — > ac/n, which tends to zero for large n, and thus the series would be convergent. 
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